Lack of pollinators selects for increased selfing, restricted gene flow and resource allocation in the rare Mediterranean sage Salvia brachyodon

Range contraction and habitat fragmentation can cause biodiversity loss by creating conditions that directly or indirectly affect the survival of plant populations. Fragmented habitats can alter pollinator guilds and impact their behavior, which may result in pollen/pollinator limitation and selection for increased selfing as a mechanism for reproductive assurance. We used Salvia brachyodon, a narrowly distributed and endangered sage from eastern Adriatic, to test the consequences of range contraction and habitat fragmentation. Molecular data indicate a severe and relatively recent species range reduction. While one population is reproductively almost completely isolated, moderate gene flow has been detected between the remaining two populations. The high pollen-to-ovule ratio and the results of controlled hand pollination indicate that S. brachyodon has a mixed mating system. Quantitative and qualitative differences in the community and behaviour of flower visitors resulted in limited pollination services in one population where no effective pollinator other than pollen and nectar robbers were observed. In this population, self-pollination predominated over cross-pollination. Various environmental factors, in which plant-pollinator interactions play a pivotal role, have likely created selection pressures that have led to genetic and phenotypic differentiation and different resource allocation strategies among populations.

to scarce pollen supply 22 .However, pollen limitation can also select for increased autogamy as a mechanism of reproductive assurance, when it increases average seed production beyond what is produced in the absence of autogamy 23 .Under these circumstances, autogamy, at least in some years, may help explain the maintenance of mixed mating systems in a plant population [24][25][26] .
It is widely accepted that generalists are less affected by habitat fragmentation than specialists and that habitat fragmentation often leads to the extirpation of specialists followed by an influx of generalist species [27][28][29] .However, the specialization of plant-pollinator interactions is asymmetric and the plant specialization cannot be considered in isolation from the degree of specialization of the mutualist partners 30 .Nevertheless, disruption of pollination and seed dispersal mutualisms involving specialists and/or keystone species is particularly problematic because of the potential cascading effect on the rest of the community 10 .One species that is characteristic of specific sub-Mediterranean rocky grassland communities in the eastern Adriatic is Salvia brachyodon Vandas (Fig. 1a).It is a narrowly distributed and endangered Mediterranean sage (Fig. 2b) that occupies dry (sub)Mediterranean grasslands (Fig. 2c and d) with highly fragmented populations 31 .We selected all currently known populations of S. brachyodon to test the hypothesis that the species suffers from habitat contraction and fragmentation.Specifically, we asked: (a) what is the extent of gene flow among populations?(b) Are there differences in floral biology, pollination ecology, and mating systems between populations?(c) Is there evidence of pollen/pollinator limitation?(d) What is the difference in degree of ecological specialization between populations?(e) Does selfpollination, as opposed to cross-pollination, lead to trade-offs in flower size and seed weight?To assess genetic structuring and gene flow among populations, we used plastom sequences and DNA fingerprinting data in all extant populations, while two of the three extant populations were used to study the mating system by in situ controlled hand pollination.We examined floral morphology and biology, seed weight, and visitor/pollinator assemblages to address the problem of resource allocation triggered by the shift from outcrossing to autogamous selfing and the extent of phenotypic differentiation and ecological specialization.Since S. brachyodon is an endangered species, our study additionally aimed to identify the source population for its possible translocation activities at potentially suitable sites.

Population genetic structuring and gene flow
From 265 analysed samples, 78 multi-locus genotypes (MLGs) were recognized.In population OR, only 11 MLGs were identified, while 36 and 31 MLGs were detected in populations PE and KO, respectively.Following the number of detected MLGs, slightly higher levels of observed and expected heterozygosity and allelic richness were detected in populations PE and KO if compared to population OR which was seemingly characterized by depleted levels of genetic variability (Table 1).Nevertheless, no departures from Hardy-Weinberg equilibrium were detected in any of the studied populations.The fixation index (F st ) values were significant for all pairs of populations.The strongest differentiation was detected between populations PE and OR (F st = 0.22), the weakest between populations KO and OR (0.14), while the fixation index value of 0.17 was detected for the populations PE and KO (Fig. 2a).Accordingly, the highest values of gene flow (N m ) were obtained between populations OR and KO (1.55), and the lowest between populations OR and PE (0.88), while value between populations KO and OR was 1. 18.The results of F st statistics and gene flow followed relative geographic distances between pairs of populations.Results of both PCA (Fig. 2b) and Bayesian assignment test (Fig. 2c) support the presence of two gene pools (genetic clusters).The highest ΔK was obtained at K = 2 (2195.7),while the ΔK at K = 3 was 203.1.The larger genetic cluster is comprised of MLGs from populations KO and OR, while genotypes from population PE form the other one.The lengths of the rpl32-trnL and rps16-trnK sequences were 796 and 685 bp, respectively.No nucleotide substitutions and indel events were detected.

Flower life-span and sexual functioning
Mean flower life-span of S. brachyodon was one day (1.15 ± 1.27, 97; mean ± SE, n) in unbagged control flowers, and up to three days in the bagged test flowers (Fig. 3a), and did not differ between populations PE (n = 45) and OR (n = 52; χ 2 = 0.81, p = 0.36).
At population PE, 41 flowers (91%) had a hole in a corolla tube, 9 already in the bud stage, 4 of which subsequently wilted immediately.The number of flowers damaged by nectar robbers was significantly lower in population OR (58%, χ 2 = 13.29,df. 1, p < 0.001), where no holes were detected in the bud stage (Fig. 4).The anthers dehisced immediately after flower opening, and pollen was generally available in the anthers until the second day, when it was almost completely depleted by floral visitors in most observed flowers.The highest pollen germination was observed on the first day when flowers opened (mean = 87%), and decreased steadily (83%) until the second day of anthesis and significantly (68%; p < 0.05) on the third day.Low stigma receptivity was observed on the first day (12%), but increased significantly (p < 0.05) with the onset of stigma bifurcation on the second (80%) and third day (85%).The data reveal weak protandry with an overlap of both functions on days two and three of flower life-span, coinciding with high pollen germinability and stigma receptivity.

Mating system, pollen limitation and pollen-to-ovule ratio
Pollination treatments significantly affected fruit set in both studied populations (χ 2 = 180.47,df.= 10, p < 0.0001; Fig. 3c, Supplementary Table S1).Spontaneous autogamy (A s ) at populations PE (8%) and OR (6%) resulted in significantly fewer fruits than any other treatment except the control treatment (C) at population PE (4%).Compared to A s , the induced autogamous treatment (A i ; performed only in population OR) produced more fruits (42%, p < 0.05).Fruit set in geitonogamous treatments (G) were not significantly different from treatments www.nature.com/scientificreports/A i and pollen limitation (PL: 64% and 69% for populations PE and OR, respectively) in both populations, but was significantly higher than treatment C in population PE and lower than the xenogamous treatment (Xe) in population OR (88%).The highest fruit set (88%) was obtained with treatments Xe and C in population OR.
With the exception of treatment A i (OR), plants in population OR had higher or significantly higher fruit set compared to plants in population PE, regardless of the treatment applied.Inbreeding coefficients (δ i ) were -0.07 in populations PE and 0.39 in population OR.The number of pollen grains produced per flower in population

Flower visitors
In total, we observed 71 individuals of 24 different insect species on the flowers of Salvia brachyodon of populations PE and OR (Table 2).Most visitors belonged to solitary (18 species) or social bees (5 species), and only one butterfly (Sphingidae) and one wasp (Vespidae, Eumeninae) species were observed.According to their behaviour, pollen thieves dominated, represented in most cases by small solitary bees (14/24 species), followed by both nectar and pollen thieves and nectar thieves (4 and 3 species, respectively).Primary nectar robbers, which gain  S3).In population PE, an individual of Episyrphus balteatus, a marmalade hoverfly (Syrphidae), was recorded with its parasite foraging for pollen between the observation periods (and thus the record not included in the analysis; Video S4).In both populations, the frequency of flower visits was highest between 10 00 and 12 00 and higher in the morning (7 00 -9 00 ) than in the evening (16 00 -18 00 ; Fig. 3d).While the number of different insect taxa visiting the flowers of S. brachyodon was nearly identical at each locality (14 and 13 in populations PE and OR, respectively), the insect fauna differed significantly in composition (χ 2 = 55.313,df.23, p < 0.001), flower visitation frequency, Shannon diversity (p = 0.02), dominance (p = 0.0008) and evenness (p = 0.01).Although S. brachyodon flowers were visited by insects in the population OR almost twice as often as in the population PE, most flowers in the population OR were visited by individuals of B. argillaceus, as evidenced by the significantly higher dominance and Shannon diversity index and the significantly lower evenness index for the population OR.Between the two populations, the flowers of S. brachyodon shared visitors (nectar robbers, thieves, and pollen thieves) from only three taxa (12.5% of total pollinator assemblage), and the only pollinator, B. argillaceus, was observed only on flowers of population OR.

Floral morphology and seed weight
Flowers between populations differed significantly in 11 of the 16 measured characters (Table 3, Fig. 5a).In general, the largest character values were observed in the OR population, while the smallest were observed in  7), greater spacing between the inner hairs of the corolla tube and the nectaries (8), and the labium and lower thecae (11).Also, the distances between the lower and upper thecae (15) were the largest in the OR population.On the other hand, the length (2) and width (3) of the sepals, the length of the corolla tube (4), and the distance between the hairs of the inner corolla tube and the nectaries (8) were the smallest in the population  PE.Most of the flower characteristics measured in the KO population were between the values obtained from the flowers in the OR and PE populations.The least variable characters were corolla length (1, coefficient of variation = 4.4-14.5%),upper theca length (14, 8.6-11.6%),corolla tube length and width (4, 8.6-18.2%and 5, 8.8-15.9%,respectively), and calyx length (2, 10.2-12.1%),while the most variable characters were the distance between the upper theca and the stigma (16, 38.2-58.4%), the distance between the stigma tip and the labium (12, 19.7-43.3%),and the length and width of the labium (9, 19.3-36.6% and 10, 14.9-25.1%,respectively).Mean values of seed weight among the different populations ranged from 2.7 to 7.4 mg (Me = 3.9 mg; min.indiv.value 0.1 mg, max.9.7 mg; Supplementary Table S2).Seed mass differed significantly between pollination treatments at both sites (χ 2 = 65.52,p < 0.001; Fig. 5b).In populations PE and OR, selfing (treatments A and G) as opposed to outcrossing (Xe) resulted in lighter seeds, although significantly higher values (p < 0.05) were obtained only in treatment Xe in population OR.Seeds from treatment G in population PE were significantly lighter (p < 0.05) than those from treatments Xe in populations PE and OR.

Discussion
The absence of polymorphisms in analysed plastome regions indicate a relatively recent and severe bottleneck event.Such a result was somewhat unexpected because in closely related S. officinalis, cpDNA loci were characterized by high polymorphism levels 32 .The ecologically similar and co-occurring Edraianthus tenuifolius Table 1 in 31 experienced range contraction during the Last Glacial Maximum, but subsequently expanded its range north-westward along the climatically buffered eastern Adriatic coast 33,34 .However, the range of S. brachyodon, which was possibly more widespread in the past, has remained restricted and without evidence of recent expansion, contributing to the fate of Mediterranean flora with extremely high rates of narrow endemism 35,36 .Indeed, historical literature data report the occurrence of the plant approximately 80 km north of population PE 37 .However, recent confirmations of this occurrence are lacking.It remains intriguing that, despite the presence of suitable habitats along the eastern Adriatic coast, albeit fragmented, the species has remained restricted to only three extant populations.
The results of the microsatellite data analyses (Fig. 2, Table 1) indicate moderate gene flow between populations OR and KO, suggesting recent distribution range and habitat fragmentation and/or the existence of some undetected bridging populations enabling gene exchange between them.At the same time, population PE remains substantially more reproductively isolated.Such a finding was additionally supported by both Bayesian assignment test and PCA results (Fig. 2b and c) that indicate the existence of two well-defined genetic clusters, one comprising populations KO and OR, and another only population PE.The continued encroachment of deciduous woody species into calcareous/dolomitic rocky grasslands and habitat loss present a threat to all S. brachyodon populations, but is particularly evident in population KO 31 , where the number of plants in sample plots for genetic analyses declined from 70 + to only a few in just five years (personal observation).An excess of homozygotes is usually a consequence of increased levels of breeding among closely related individuals.However, our data suggest that fragmented and spatially isolated population KO did not suffer from reduced genetic diversity or an excess of homozygotes.It appears that even moderate pollen-mediated gene flow can contribute to long-term population viability, as demonstrated in Ziziphus lotus (Rhamnaceae), a sclerophyllous shrub that occurs in semiarid habitats throughout the Mediterranean region 38 .Nevertheless, habitat fragmentation and its subsequent loss, followed by a decline in population size that affects mate availability via loss of genetic diversity, will likely negatively impact the long-term persistence of S. brachyodon 39 , necessitating active conservation measures especially for population KO.In such cases, population size, genetic diversity, and geographic distance have often been used to make decisions.However, Lawrence and Kaye 40 suggested that ecological similarity is a better predictor of plant performance and establishment when identifying source populations for reinforcement and translocation actions.Community assemblage analysis of all extant populations showed the greatest similarity between populations OR and KO, while population PE differed the most 31 .Combined with the lower degree of differentiation, caused by moderate gene flow between populations OR and KO found in the present study, population OR appears to be a suitable source population for reinforcement and translocation actions in population KO and neighbouring suitable habitats.
Flower longevity did not differ between populations of S. brachyodon, despite considerable differences in site ecology and nectar robber's abundancy and activity.Similar to S. gesnerifolia Lindl.& Paxton 41 , nectar robbers were much more active in drilling holes in young flowers, especially in the population PE, but did not affect flower longevity (Fig. 4).Consistent with the minimum longevity hypothesis 42 that pollination-induced senescence should be more common in protandrous species 43,44 , unmanipulated flowers (treatment A s ) of S. brachyodon that were protected from insect visitors (Fig. 3a) lasted much longer than open-pollinated (C and PL) or handmanipulated flowers (A i , G and Xe).A similar observation was made by Aximoff and Freitas 45 in S. sellowiana Benth., where bagged and non-manipulated flowers remained functional twice as long as manipulated flowers.
During anthesis, S. brachyodon produced nectar during the day and night in a pattern, quantity, and concentration very similar to that of S. pratensis 46 .In general, as temperature increases, the rate of nectar secretion also increases 47 .However, in the afternoon, the amount of nectar secreted by flowers of S. brachyodon decreased, probably as a result of heat and water stress 48 .This pattern was more pronounced at population PE, where the decrease in nectar secretion was observed as early as midday (Fig. 3b), when plants appeared wilted but recovered during the night.
Significant differences in the composition and abundance of flower visitors between the two sites (Table 2) were not surprising.However, temporal and spatial variations in the quality and quantity of flower visitors (including pollinators) appear to be the rule rather than the exception [49][50][51][52] .Despite the fact that results of different studies on the composition and abundance of insect flower visitors in Salvia are difficult to compare due to different sampling methods and duration of observation periods 53 www.nature.com/scientificreports/large number of insect taxa, of which only B. argillaceus effectively triggers the lever mechanism, which, despite the wide variation in flower visitor composition, suggests that S. brachyodon is an ecological specialist in which pollination niches are defined by a major pollinator.Consequently, any change in this flower visitor is usually accompanied by a change in the pollination system 59 .Changes in the pollination niches of specialized plants associated with replacement of the main pollinator usually result in adaptive pollinator-mediated divergence in floral traits [60][61][62] .The fruit set observed in the field and the results of the hand pollination experiment indicate that in S. brachyodon, although closely adapted to a single pollinator species, thieves and robbers of nectar and pollen occasionally contribute to pollen transfer, allowing sexual reproduction in the absence of the pollinator.Of the various life history traits (including the number of fruits and seeds produced), mating systems have the greatest influence on patterns of polymorphism 63 .The high pollen-to-ovule ratio and the results of controlled hand pollination indicate that S. brachyodon has a mixed mating system and generally prefers cross-pollination to self-pollination (Fig. 3c).When pollinators were excluded (A s ), fruit set declined significantly, suggesting that this species requires pollination vectors to ensure reproductive success.Fruit set after autogamous (A i ) and allogamous treatments (G, Xe, and PL) did not significantly in either population, except for the xenogamous treatment (Xe) in population OR, indicating that individuals in this population clearly prefer outcrossing to selfing.Plant mating systems may vary considerably among taxa and populations both spatially and temporally [64][65][66] due to a variety of intrinsic and extrinsic factors that often have confounding effects.The low fruit set in control treatment in PE is best explained by a lack of pollinators (Table 2), but can also be a result of nectar robbing (Videos S1 & S2) and herbivory.The effects of nectar robbers are very complex and can affect plant reproductive success in several ways 67,68 .Our results suggest that nectar robbers may have a direct or indirect detrimental or neutral effect on fruit set in S. brachyodon.In population PE, four out of nine flower buds (44%) with boreholes drilled by nectar robbers wilted before anthesis, which alone could reduce fruit set by 8% (Fig. 4).In the same population, severe habitat destruction caused by wildfire may have resulted in significant turnover of flower visitors.Here, nectar robbers may have disrupted the pollinator niche that hasn't recovered after the fire, resulting in reduced fruit set, as shown in the ecologically similar S. fruticosa Mill. 54,69.Sparse fruit set in population PE observed in the control treatment and in neighbouring plants not included in the experimental design suggests that thieves and robbers may also contribute directly to pollination success, albeit to a negligible extent.
Flower characteristics differed significantly among populations (Fig. 5a, Table 3), which could be partly explained by differences in site ecology and composition of flower visitors.On the other hand, mating systems have the greatest influence on patterns of polymorphism.Increased isolation between populations results directly from selfing and indirectly from evolutionary changes, i.e., smaller flowers 63 .According to this scenario, selfing populations are more genetically differentiated than outcrossing populations and reduce allocation to attraction 70 .In all controlled pollination treatments, except PL in population OR, seed weight was higher in population OR than in genetically differentiated population PE, although significantly only in treatment Xe in population OR (Fig. 5b).These results are also consistent with the fruit set (Fig. 3c).In general, allogamous treatments in S. brachyodon resulted in higher seed weight than autogamous treatments, which was also demonstrated in S. verbenaca 71 .On the other hand, seed weight and fruit set in treatments A and G in population PE were not significantly different from treatments PL and C (Figs. 3d and 4), suggesting that most of the pollen flow in population PE occurs between flowers of the same plant.Morphometric analysis of the two populations revealed significantly higher values for vegetative traits (plant height, inflorescence length, number and length of internodes, number and length of primary branches) and lower values for reproductive traits (calyx tube and teeth lenght) for population PE in comparisson to population OR 72 .This is consistent with the results of our morphometric analysis of flowers, where the higher values for flower traits were obtained in population OR, suggesting different strategies of resource allocation between vegetative and reproductive parts in the two populations.

Conclusions
The results of this multidisciplinary study, which included different types of genetic markers and anthecological analyses, suggest that Salvia brachyodon is an ecologically specialised sage exposed to persistent and extensive environmental pressure.Range contraction and habitat fragmentation during the evolutionary history of the species have had a significant impact on present-day genetic structuring, which is associated with strong differences between populations in reproductive biology and survival strategies.Results derived from plastom DNA analysis indicate a recent and severe genetic bottleneck that has erased the majority of the species' genetic diversity.The extant populations have found themselves in habitats that are isolated geographically, reproductively, and genetically.Various ecological factors, in which plant-pollinator interactions play a pivotal role, have likely created selection pressures that have led to genetic and phenotypic differentiation.This recent differentiation, detected by microsatellite markers and coupled with contrasting characteristics in the pollination biology and resource allocation strategies of the species, sheds light on the ongoing and dynamic process of its adaptive divergence.

Study system
Among the European sage species, Salvia brachyodon (Fig. 1a) stands out for its large flowers, along with S. eichlerana Heldr.ex Halácsy and S. ringens Sibth.& Sm.However, it is noteworthy that S. brachyodon has relatively shorter calyx teeth compared to the other two species.Morphologically, S. brachyodon closely resembles S. ringens, an allopatric endemic of the southeastern Balkans.Phylogenetic studies have revealed that both S. brachyodon and co-occurring S. officinalis L. belong to the Salvia subgenus Salvia 73,74 .Compared to the evolutionarily closely related and sympatric S. officinalis, S. brachyodon prefers relatively cooler, moister, deeper, and more nutrient-rich soils that occur on dolomite or dolomitic limestone at higher elevated sites.Their flowering periods also do not overlap.Nowadays, there are only three known populations of S. brachyodon, all located in www.nature.com/scientificreports/ the central part of the coastal Dinaric Alps along the eastern Adriatic Sea (Fig. 1b), all with a small number of individuals: Pelješac Peninsula, Konavle and Mt.Orjen 31 , hereafter referred to as populations PE, KO and OR, respectively.Despite its ability to cope better than S. officinalis with interspecific competition through clonal reproduction 75 , populations of S. brachyodon are severely threatened by the abandonment of traditional land use and by fire, making the species endangered (EN) according to IUCN criteria 31 .In fact, population KO is on the brink of extinction (personal observation).Although specimens of S. brachyodon appear fairly uniform morphologically, Liber et al. 72 found significant differences in some vegetative and reproductive traits among populations.Individuals in population PE are generally more robust, have a greater number of primary branches and internodes, and longer inflorescences, whereas individuals in population OR have longer calyces.From the point of view of population genetics, population PE was studied in detail by conducting a high-resolution spatial genetic analysis, confirming a high degree of clonality within the population 75 .In 1998, this population suffered a severe wildfire that completely changed the habitat characteristics and community assemblage.Contrary to expectations, a study conducted 15 years after the fire at the same site, showed a slight heterozygote excess, a high level of genetic diversity, as well as clonal reproduction and a genetic bottleneck 76 .Nowadays, S. brachyodon, with a peak flowering period in late July and mid-August, is a keystone species and an edificator in four floristically and ecologically well-defined groups of stands distributed in three localities.Although these stands represent different syntaxa, they all belong to the dry eastern (sub)Mediterranean rocky grasslands, which were once more widespread and continuously distributed (Fig. 1c and d).However, their current state varies in terms of successional stages towards (sub)Mediterranean forest vegetation 31 due to abandonment of traditional landscape management.Population OR occupies stony pastures and grasslands 77 characterized by shallow calcareous soils derived from Jurassic dolomitic limestone 78,79 .In addition to these habitats, this population can also be observed in fringe vegetation and thermophytic forests containing Quercus pubescens Willd., Carpinus orientalis Mill., and Fagus sylvatica L. 80 .Similarly, population PE is not exclusively confined to stony pastures.It historically thrived on humus-rich soils within stands of Dalmatian Pine (Pinus nigra subsp.dalmatica (Vis.)Franco), which are situated just above the climazonal evergreen Mediterranean vegetation 81,82 .The geological composition of the bedrock supporting population PE consists of dolomites and limestones dating back to the lower and upper Cretaceous periods 83,84 , while population KO is primarily found on Jurassic dolomitic limestone (ibid.).The mean annual temperature in S. brachyodon sites varies, ranging from 12-13 °C in populations PE and KO to 10-11 °C in population OR 85 .The amount of annual precipitation increases in a west-east direction, with populations PE, KO, and OR receiving approximately 1300-1400 mm, 1500-1750 mm 86 , and more than 3700 mm 87 , respectively.The precipitation regime in these areas follows a Mediterranean pattern.

Population genetics
In total, leaf material from 265 individuals was collected in 2015: 119, 70 and 76 from populations PE, KO and OR, respectively.After rapid desiccation of the plant material in silica gel, DNA was extracted using GenElute plant genomic DNA miniprep kit (Sigma-Aldrich, St. Louis, MO, USA).The set of analysed loci, the PCR procedure, and genotyping, were as described in Radosavljević et al. 75 .To identify multi-locus genotypes (MLG) and to assess the levels of clonality (genotypic richness (R)), RClone package version 1.0.2 88was used.After the identification of individuals (i.e., ramets) that shared the same genotypes and before the further analysis of populations' structure and diversity, replicates were removed from the dataset.After the identification of the MLGs, to assess the phylogeographic structure of the species, two regions of the plastome genome were analysed, rps16-trnK and rpl32-trnL.These regions were selected based on recently published results 32 , where they were characterized by significant variability in the closely related S. officinalis.For this purpose, 20 MLGs from populations PE and KO, and 10 MLGs from population OR were randomly selected.PCR reactions were performed in a total volume of 15 μL containing 8 μL Nuclease-Free Water (Qiagen®), 1.5 μL 10 × PCR buffer (TaKaRa Taq™ Hot Start Version), 1.2 μL dNTP solution (2.5 mM each, TaKaRa TaqTM Hot Start Version), 0.6 μL 5 μM of each forward and reverse primer, 0.1 μL TaqHS polymerase (5 U/μL, TaKaRa TaqTM Hot Start Version) and 3 μL of 1 ng/μL DNA.The PCR conditions were as follows: 94 °C for 5 min, followed by 30 cycles of 95 °C for 1 min, 50 °C for 1 min and 65 °C for 4 min, and a final elongation step of 5 min at 65 °C.Exonuclease I (Thermo Scientific™) and FastAP Thermosensitive Alkaline Phosphatase (Thermo Scientific™) was used for the purification of PCR products.The sequencing of the amplicons in both directions was performed using an ABI 3730XL analyser (Applied Biosystems).Sequences were cleaned with Codon Code Aligner v. 9.0.1 (CodonCode Corporation, Dedham, USA) and aligned using the sequence alignment program ClustalW implemented in the MEGAX software 89,90 .Sequences were deposited in GenBank.

Flower life span and sexual functioning
Flower longevity was investigated by marking a total of 97 flowers, with 45 and 52 from populations PE and OR, respectively, on multiple individuals prior to their opening.Subsequently, the newly opened flowers were observed daily at 8 00 , 14 00 , and 19 00 h until senescence, and the duration of each flower was recorded in hours.
The presence and first appearance of the hole in the corolla tube just above the calyx teeth caused by nectar robbers, a damage that could affect flower longevity, was recorded for different stages of flower anthesis.We assessed pollen viability in upper and lower thecae throughout the flower lifespan using diaminobenzidine reaction (DAB) 91 , assuming that viable pollen, unlike dead and aborted pollen, turns dark reddish brown due to enzyme activity 92 .The effectiveness of the staining was first tested on pollen killed by high temperature.When evaluating the percentage of pollen stained, at least 200 grains (and up to 400) per flower were counted.Since the stigma viability in Salvia generally correlates with the opening of the stigma lobes 93 , we first calibrated the intensity of staining of the stigma receptive surface using DAB test and stigma bifurcation using five grade scale for staining (1-stigma not stained, 2-stained only the tips of stigma, 3-at least half of the stigma stained, 4-more than half of the stigma stained, 5-stigma fully stained) and a three-point scale (1-stigma closed, 2-stigma lobes half opened, 3-stigma lobes fully opened) for stigma opening.Stigma receptivity was then assessed by using results of both staining and opening of lobes.A total of 82 and 22 flowers were used to evaluate pollen viability and stigma receptivity, respectively.Both pollen viability and stigma receptivity were tested in population OR.

Floral rewards
Based on our initial field observations, it became evident that flower visitors of S. brachyodon engage in foraging for nectar and/or pollen.Consequently, we proceeded to quantify both aspects.To assess daily nectar production, we selected a random sample of flower buds (69 from population PE and 37 from population OR) belonging to at least 20 different genets per population.This selection aimed to minimize potential damage caused by insect visitors.The chosen flower buds were individually enclosed in tulle bags 24 h before their expected opening, ensuring the exclusion of nectar collection by insect visitors.Subsequently, the bagged flowers were carefully opened, and nectar was extracted using 5 µl microcapillaries (Drummond) at three different time intervals: 6 00 -7 00 h, 12 00 -14 00 h, and 17 00 -19 00 h.The length of the nectar column was measured using digital callipers.
Once the nectar collection was completed, the bagged flowers were resealed to prevent further nectar removal or contamination.To determine nectar concentration, the collected nectar was dispensed onto a handheld portable refractometer (0-50°Brix, Eclipse Professional, Bellingham and Stanley Ltd), and the % Brix (i.e.g sucrose per 100 g solution) was recorded.The nectar standing crop was determined by measuring a total of 31 well-preserved and freshly opened flowers (from the population PE) at 8 00 h and 19 flowers (from the population OR) at 14 00 h.These flowers were randomly selected from several genets per population to ensure representativeness.

Mating system, pollen limitation and pollen-to-ovule ratio
Pollination ecology and floral biology studies on Salvia brachyodon were performed in situ only in populations PE and OR due to a low number of individuals in population KO.Buds selected for treatments were covered with a tulle to prevent insect interactions.The following treatments were applied to determine the mating system, effects of insect exclusion and pollen source on fruit production: (a) spontaneous autogamy (A s )-individual flowers covered with a tulle were left intact; (b) induced autogamy (A i )-individual flowers were pollinated with their own pollen, then emasculated by cutting of the anthers and covered again with the tulle; (c) geitonogamy (G)-individual flowers were carefully emasculated in the bud, pollinated with pollen from flowers of the same plant and then covered again with a tulle; (d) xenogamy (Xe)-individual flowers were emasculated in the bud, pollinated with a fresh pollen mixture collected from several different plants from the same population and covered again with a tulle; (e) supplementary pollination (PL)-flowers were pollinated with a fresh pollen mixture collected from several individual plants and left uncovered; (f) control (C)-flowers were labelled and left intact.A total of 288 (PE) and 175 (OR) flowers were included in the treatments.After 30 days, fruit production was recorded.From the results of hand pollinations, we calculated index of inbreeding depression (δ i ) related to the mating system.To distinguish possible effects of dichogamy from genetic incompatibility, autogamous performance was determined from the results of geitonogamous pollen transfer.Hence, the degree of inbreeding depression was determined by the ratio between fruit set and reproductive success of geitonogamously (wG) and xenogamously pollinated flowers (wXe), as suggested by Charlesworth and Charlesworth 94 : δ i = 1 − (wG/wXe).
To estimate the pollen-to-ovule ratio (P/O ratio) and overall pollen production, we counted pollen using a Fuchs-Rosenthal counting chamber; we took 10 µl of the pollen suspension from a single thecae and filled the chamber with a volume of 0.2 µl.Pollen grains were counted and then the number of pollen grains per anther and flower was calculated.We performed 10 replicates with one upper and one lower thecae per flower.P/O ratio as the total number of pollen grains produced per flower divided by the number of ovules was then calculated as a proxy for the breeding system 95 .The specific number of manipulations in studies on sexual functioning, floral rewards and mating system is indicated in Figs. 3 and 5.

Flower visitors
The assemblage and behaviour of flower visitors were assessed by direct observations, video recordings, and net catches during the whole time of investigation from July 20 to July 27 2016 in population PE and from July13 to July 19 2016 in population OR.If necessary, visitors were captured for identification and subsequently preserved in ethyl acetate for later dissection and identification in the laboratory.Voucher specimens are deposited at the Natural History Museum Rijeka.Observation periods were determined based on the previous years' experience of insect behaviour, visitation frequency, and dynamics during the day.Quantitative observations, involving the count of visitors, as well as qualitative observations, documenting the specific taxa visited and their corresponding behaviors, were conducted within 4 m 2 quadrats with homogenous plant cover.These observations took place on sunny and windless days during the main flowering period, specifically between 7 00 -9 00 , 10 00 -12 00 , 13 00 -15 00 and 16 00 -18 00 h GMT.The observations were carried out on July 13-14, 2017, in the OR population and on July 20-21, 2017, in the PE population, resulting in a total of 32 h of observation.Flower visitors were classified as legitimate visitors if they entered the corolla through the corolla tube, or illegitimate if they accessed pollen or nectar by other means.Visitors were further classified as pollen (PT) or nectar thieves (NT) if they legitimately entered the flower to search for pollen or nectar without touching the anthers and/or stigma, and as nectar robbers (NR) if they gained access to nectar through a hole in the corolla tube 96 .Legitimate visitors who effectively triggered the staminal lever mechanism and contacted both the thecae and the stigma surface were classified as pollinators (P).www.nature.com/scientificreports/

Flower morphology and seed weight
To examine potential resource allocation among reproductive parts between populations, we measured 16 quantitative floral traits in situ using a digital calliper (see Table 3) to estimate flower size and functioning depending on morphology.The selection of specific traits for measurement can have a significant impact on pollinator composition and efficiency in Salvia spp. 56,58.Seeds obtained from controlled hand pollination in populations PE and OR were weighed to the nearest 10 μg (Ohaus AP250D) in order to check whether seed weight depend on pollen source (selfing vs. outcrossing) and for the differences between populations.Voucher specimens are deposited in the herbarium of the Natural History Museum Rijeka (NHMR).

Statistical analyses
Population genetic datasets were analysed using different statistical tools.To estimate observed and expected heterozygosity and inbreeding coefficient, we used GENEPOP software 97 .R packages "PopGenReport" 98 and "hierfstat" 99 were used to assess the allelic richness and pairwise F st values, respectively.For the assessment of the genetic structure, in order to elucidate the structuring patterns and to make certain assumptions about the data, the principal component analysis (PCA) was performed as implemented in the R package "adegenet" 100 , followed by the Bayesian assignment analysis using STRU CTU RE software ver.2.3.3 101 .Thirty runs per cluster (K), with K ranging from 1 to 4, were performed.Each run consisted of a burn-in period of 200 000 steps followed by 1 000 000 MCMC replicates.To process the obtained results and to select the most likely number of clusters following Evanno et al. 102 , STRU CTU RE HARVESTER v0.6.92. 103was used.Finally, runs were clustered and averaged using CLUMPAK 104 .Gene flow between populations (N m ) was calculated from the pairwise F st values using the formula: N m = (1 − F st )/4 F st 105 .The data on floral biology and mating system were subjected to analysis using generalized linear models (GLM) with fruit set as the response variable and treatments and populations as predictors, considering all populations and each population separately.The response variable was fitted to a binomial distribution with a logit link function.Likelihood ratio test was performed to compare the full model with a restricted model and by calculating p values using the χ 2 distribution.Differences between levels of each effect were analysed post hoc by multiple comparisons of means with Tukey contrasts, adjusting data for normality and testing for homogeneity of variance.All statistical analyses were performed using the R version 3.1.1 106and packages 'stats' , 'car' , and 'multcomp.
Differences in stigmatic receptivity and pollen germinability along the flower lifespan were investigated using GLM with the response variables (staining/no staining for stigma receptivity and pollen germinability) being adjusted to a binomial distribution and a logit link function.
Nectar volume and sugar concentration were adjusted to a Gaussian distribution with an identity link function.
Data on mating system were analysed with fruit set as the response variable and manipulations as predictors, considering all populations and each population separately.The response variable was fitted to a binomial distribution with a logit link function.
We examined the dominance, evenness, and diversity of pollinators across populations, employing a range of indices.To compare these indices, we employed a permutation test for equality.Additionally, we assessed differences in flower visitor fauna using the χ2 test implemented in PAST 107 .Results were considered significant if the probability for the null hypothesis was below 0.05.
To analyse the differences in flower and seed characteristics, we performed the Kruskal-Wallis test, as the data sets did not meet the assumptions of normality and homogeneity of variance.This non-parametric test allowed us to assess the equality of medians.Subsequent post-hoc pairwise comparisons were conducted through the Mann-Whitney test using PAST.Statistical significance was determined if the probability for the null hypothesis fell below 0.05, with p-values adjusted using the Bonferroni correction for multiple comparisons.

Ethical approval
Formal identification of the plant material (Salvia brachyodon) was provided by Boštjan Surina.Voucher specimens of the plant material used in the study have been deposited in a publically available herbarium at the Natural History Museum Rijeka (Index Herbariorum acronym NHMR: https:// sweet gum.nybg.org/ scien ce/ ih/ herba rium-detai ls/? irn= 157889).Deposition (voucher) numbers are included in the manuscript (Table 1).Permission to perform the research and to collect voucher specimens of Salvia brachyodon was issued by Nature Protection Directorate of the Ministry of Economy and Sustainable Development of the Republic of Croatia (UP/I-352-04/22-08/15, 517-10-1-1-22-4).All the experimental research and field studies on study plant complied with relevant institutional, national, and international guidelines and legislation.

Figure 2 .
Figure 2. Population genetics of Salvia brachyodon.(a) Approximate geographic positions of studied populations of Salvia brachyodon: Pelješac peninsula-PE, Konavle-KO, Mt.Orjen-OR, with indicated distance (km) and gene flow (N m ) between populations.(b) Principal component analysis (PCA) using the first and the second axis of specimens of all extant populations based on microsatellite data.(c) Genetic structure and assignment of individuals into classes as assessed by the computer program Structure.Each individual specimen is represented by a single vertical line; each colour represents a cluster, and the length of the coloured segments indicates the individual's estimated proportions of membership in those clusters.Population's symbols (enlarged) in (a) follow those in (b), while arrow thickness in (a) indicate intensity of gene flow.Colour codes in (b) are congruent with the results of genetic structuring (c).

Figure 3 .
Figure 3. Flower biology, mating system and insect visitation of flowers of Salvia brachyodon.(a) Sexual functioning of flowers of Salvia brachyodon from population Pelješac peninsula-PE.(b) Nectar volume and concentration produced by flowers during the day in population PE and population Mt.Orjen-OR.(c) Controlled hand pollination experiments conducted to study the mating system and to assess pollen limitation for populations PE and OR.Pollination treatments: A s -spontaneous selfing, A i -induced selfing, G-geitonogamy, Xe-xenogamy, PL-pollen limitation, C-control.(d) Mean insect visitation frequency of flowers during the 2 days in four observation periods in populations PE and OR.Different letters indicate statistically significant differences in male and female function according to flower age (a) and nectar volume and concentration (b; p < 0.05).Numbers in parentheses indicate the number of flowers manipulated (a-c).